Effect of material stiffness on hardness: a computational study based on model 
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We investigate the dependence of the hardness of materials on tlieir elastic stiffness. This is 
possible by constructing a series of model potentials of the Morse type; starting on modelling natural 
Cu, the model potential exhibit an increased elastic modulus, while keeping all other potential 
parameters (lattice constant, bond energy) unchanged. Using molecular-dynamics simulation, we 
perform nanoindentation experiments on these model crystals. We find that the crystal hardness 
bJQ, scales with the elastic stiffness. Also the load drop, which is experienced when plasticity sets in, 

^ ' increases in proportion to the elastic stiffness, while the yield point, i.e., the indentation at which 

I plasticity sets in, is independent of the elastic stiffness. 
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^ ■ I. INTRODUCTION 

While the elastic properties of solids arc weU understood, and their description in terms of interatomic potentials is 
well established, plastic deformation of solids - and the concurrent phenomena of dislocation generation and motion 
- are more complex and their modelling presents greater difficulties. One reason certainly is that plasticity involves 
strongly non-equilibrium states in the solid, in which the material is stressed so far that interatomic bonds are broken 
and new bonds are formed. 

', The simplest quantitative measure of plasticity is given by the material hardness, i.e., the pressure with which the 
^ ■ material withstands plastic deformation. It has long been known that the hardness of a defect-free, ideal crystal, 
Q I the theoretical strength, is proportional to the material's shear modulus. This important result, which dates back to 
O . Frenkel,^'^ has been derived by considering the shear stress necessary to induce slip in a perfect lattice; modern ab 
initio quantum mechanical calculations have confirmed this result. However, the question remains how the elastic 
stiffness of a material influences dislocation generation and the emergence of plasticity in a more complex and realistic 
^ ■ situation. 

On , We choose a nanoindentation scenario to investigate the onset of plasticity, and its dependence on the material 
' stiffness. A molecular-dynamics simulation allows to provide atomistic insight into the reaction of a material to an 
] \ applied load, to calculate the force-depth curve and to extract the contact pressure and the material hardness. It 
. 1 furthermore allows to describe in detail the induced damaged patterns in the form of dislocation loops and stacking 

00 ' faults formed. 

, Our study is based on the Morse interatomic potential. While it is well known that metals cannot be described 
• in all details by pair potentials,^ this class of potentials readily allows to generate a series of potentials, which all 
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describe the same material, in which, however, exactly one property is arbitrarily changed. Since the prime aim of 
our study is to inquire into the generic dependence of plasticity and hardness on the elastic stiffness of the material, 
rather than to describe one particular material as accurately as possible, our choice of a Morse potential appears 
appropriate. We note that the Morse potential has been used previously to describe dislocations in metals. 

II. METHOD 
A. Potentials 

We use the Morse potential 

V{r) = D {exp [-2a(r - vq)] - 2 exp [~a{r ~ ro)]} (1) 

to model our material. It is characterized by three parameters: the bond strength D, the equilibrium bond distance 
ro, and the potential fall-off a. As is well known, ^'^"'^ the Morse potential cannot describe all characteristics of bonding 
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in metals; however, as discussed elsewhere, it can give a reasonable description of the elastic^'^'"'^ and also plastic 
processes occurring under indentation. We adopted this potential, since it allows to fit in a transparent way the 
materials properties. Since the Morse potential contains three parameters, it is possible to fit it to three materials 
properties; these are traditionally chosen as the lattice constant a, the cohesive energy -Ecoh, and the bulk modulus 
B. 

For example, Cu^^ with a = 3.615 A, Ecoh = 3.54 eV, and B = 134.4 GPa can be described by a Morse potential 
with D ~ 0.337 eV, tq = 2.89 A, and a = 1.33 A~^. We note that here and in the following, we cut off the potential 
at i^cut = 2.5a = 9.0375 A - i.e. including 248 neighbours - , and shift the potential to zero at this distance. 

We create a series of Morse potentials, in which the lattice constant and cohesive energy are kept unchanged, but 
the bulk modulus is set to a preassigned value. These potentials can hence be considered as describing a series of 
pseudo-Cu materials with identical bond strength, but changed bulk moduli. Table 1 reproduces the values of the fit 
parameters obtained for the potentials employed in this study. We use a Levenberg-Marquardt based optimization 
for parameter fitting; the values of a (i?coh) are reproduced within 0.1 % (1%). 

The linear elastic properties of an fee crystal are given by 3 elastic constants, cn, ci2, and C44. However, for a pair 
potential, it is always ci2= C44, so that only two elastic constants describe the elastic behaviour. We chose the bulk 
modulus 



B = 211±^ (2) 



and the average shear modulus 



^ _ Cii + 2C44 - C12 , , 

G - (3) 

to describe the elastic properties. For Cu, it is G = 39.8 GPa. Table 1 shows the shear moduli obtained for 
pseudo-Cu, and Fig. 1 demonstrates that G increases quite linearly with B; only for the highest moduli, G increases 
superlinearly with B. So in general, we may say that the entire elastic behaviour (the elastic stiffness) of pseudo-Cu 
changes in correspondence with B; we shall talk of weak and strong materials. 



B. Generalized stacking fault energy 

The generalized stacking fault (GSF) energy can be used to characterize the behaviour of a material with respect 
to formation of stacking faults, and hence dislocation formation and glide. ^^^^^^^ For its definition, wc consider an 
fee crystal, whose upper part has glided along a (111) plane with respect to its lower part by a definite amount; in 
the present paper we only consider glide vectors along the [112] direction. The GSF energy is the potential energy 
per surface area of the deformed crystal; when determining this energy, relaxation of the crystal vertical to the (111) 
plane, but no relaxation or reconstruction within this plane is allowed for.*"'^^ In our calculation, it proved necessary 
to choose the crystallite rather large, 10 lattice in lateral directions and 25 lattice constants in both vertical directions; 
we employed the conjugate-gradient technique to relax the crystal in vertical direction. 

Fig. 2 displays the variation of the GSF surface along the [112] (111) displacement. The minimum at displacement 
corresponds to fee stacking, here 7 = by definition. The second minimum is found for a displacement corresponding 
to a partial Burgers vector i[li2]a = a/V6 = 1.5 A. It describes a stacking fault with energy 7^. For our potentials, 
extremely small values were obtained, which varied in the range of 7^ = (—2 • • • + 2) mJ/m^. For comparison, the 
experimental value of 7s for Cu is 45 mJ/m^;^^'^° this value is retrieved by ab initio calculations.'^'^ We note that in 
previous investigations, a sufficiently large value of the stacking fault energy could only be obtained by adapting the 
value of the cut-off radius to i^cut = 2.2a;^ for larger i?cut, 7s strongly decreased. 

We found that the exact determination of 7s is nontrivial, since large crystallites have to be relaxed; sometimes only 
a local and not the global energy minimum may have been found. Furthermore, it is known that the stacking fault 
energy may depend sensitively on the cut-off radius of the potential.^ Hence we conclude that within the Morse pair 
potential approximation, the stable stacking fault energy of fee metals (of Cu, at least) is strongly underestimated. 
It does not - or only negligibly - depend on the elastic stiffness of the material. 

The energy barrier between the fee crystal position and the stable stacking fault is called the unstable stacking fault 
energy, 7„. It depends strongly on the elastic stiffness. As Fig. 3 shows, it increases roughly linearly with the bulk 
modulus B, and hence also with the shear modulus G, cf. Fig. 1. This may be understood since the displacement of 
two (111) planes in the crystal corresponds to a gliding motion; thus the barrier to gliding, 7„, should be strongly 
correlated with the shear modulus. 
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C. Indentation 



In order to model the plastic deformation of our material, we performed indentation simulations using the a 
appropriately adapted version of the LAMMPS molecular dynamics code.^^ Our target consists of an fee crystallite 
with a (100) surface; it has a depth of 25 nm and a square surface area of 621.2 nm^; it contains 1325598 atoms. 

The indenter is modelled as a soft sphere. We chose a non-atomistic representation of the indenter, since we are 
not interested in the present study in any atomistic displacement processes occurring in the indenter, but only in the 
substrate. The interaction potential between the indenter and the substrate atoms is characterized by a hard core of 
radius R surrounded by a softly repulsive potential^^ 



since we are in this generic study not interested in the complexities introduced by adghesion phenomena. The 
indenter stiffness was set to = 3 eV/A'^. Our indenter has a radius of i? = 8 nm. Indentation proceeds in the 
so-called velocity- controlled approach, in which the indenter proceeds with a fixed velocity, v = 20 m/s in our case, 
into the substrate. 



III. RESULTS 



A. Force-displacement curves 



Fig. 4 shows the basic result of the simulation, the force-displacement curves. In all these curves it is seen that the 
force F increases monotonically with the displacement d into the substrate until a depth tiyieid which is of the order 
of 9.36 ± 0.48 A and where the force suddenly drops. This corresponds to the onset of plastic deformation inside the 
material. The first, monotonically increasing part is due to elastic deformation of the substrate. Hertz^'^ calculated 
that for an elastically isotropic solid it holds 



F = ^Erd^^^VR. (5) 

In this relation, a single materials parameter, the so-called reduced modulus Er describes the materials elastic 
response. ^^'^^ For a rigid indenter, it may be expressed in terms of the Young's modulus E and the Poisson ratio ly of 
the substrate as 



Hertz also determined the contact pressure p; it is defined by the ratio of the normal force F divided by the contact 
area projected into the surface plane, A. We calculate A from the ensemble of contact atoms of the indenter, which 
are farthest away from the indentation axis, and approximating it as an ellipse. 

We test this law in Fig. 4b. To do this, we need to calculate the reduced elastic modulus for a (100) surface. 
Here, Young's modulus has to be calculated for deformation in (100) direction, perpendicular to the surface, while 
the Poisson ratio has to be averaged over all directions perpendicular to the (100) directions. The corresponding 
formulae are rather complex and arc provided in Rcf. 26. Fig. 4b shows that in the range of G = 80 — 180 GPa, 
the elastic behaviour follows quite uniquely the (generalized) Hertz law (5). Only for the extreme cases of G = 39 
and 278 GPa, the normalized forces F/Er are too small. We believe that this deviation from the Hertz law is due 
to (i) a poorer quality of the potential to describe the materials reaction (ii) in particular in the case of the weaker 
pseudo-Cu, G = 39 GPa, the softness of the substrate makes the fluctuations in the response of the substrate to the 
constant-velocity indentation more sizable. We note that when normalizing F/G, an equally satisfactory uniformity 
of the curves is achieved; we do not show this plot, since the theoretical foundation appears to be missing. 

At the yield point the force drops suddenly due to the onset of plasticity; this well-known phenomenon is called 
the load drop. In our series of simulations we saw that the exact positions of the yield point dyidd does not show a 
monotonic trend with the clastic stiffness of the substrate; rather dyidd fluctuates. We believe that in this constant- 
velocity indentation, the yield point may be subject to fluctuations in the simulation procedure. We therefore conclude 



4 



that within the Hmits of the fluctuations the position of the yield point does not depend on the elastic stiffness of the 
target. 

However, the load drop shows a clear increase with the elastic stiffness of the material. This is understandable since 
the atomistic reason for the load drop is the nucleation of a stacking fault in the material. Its generation requires an 
energy which scales with the unstable stacking fault energy, 7„, which has to be found to scale well with the elastic 
stiffness. 



B. Hardness 

The force F divided by the (projected) contact area A of the indenter defines the contact pressure p. After the 
onset of plasticity, it gives the hardness of the material. Fig. 5 shows the simulation results for the contact pressure. 
In view of the good quality of the scaling observed before with the material stiffness, we present our data scaled with 
the shear modulus G. This is also motivated by the consideration that the theoretical shear strength r of the material 
is given by 

T = eG, (7) 

where e is a constant, which depends on the crystal structure of the solid but is otherwise quite material-independent. 
Frenkel estimated e = l/(2\/27r) = 0.11 for the important (112) {111} glide system of fee solids.^ Recent more refined 
estimates based on density functional theory give a slightly reduced value, e — 0.085.'^"^ 

Tabor showed that the hardness measured as the contact pressure during nanoindentation amounts to 3r,^'^^ since 
the nanoindentation acts as a 'lens focussing the stress in a small volume beneath the contact point. Using Eq. (7), 
we hence expect 

H = 3t = 3eG. (8) 

Fig. 5 shows that after the onset of plasticity, the hardness is indeed proportional to the shear modulus and well 
described by Eq. (8) with a coefficient of 0.2 — 0.25, i.e., e = 0.07 — 0.08, in satisfactory agreement with the estimates 
of e given above.. The figure also appears to indicate that the fluctuations in the fully plastically regime are smaller 
for soft materials, which appears plausible. 

Note that also in the elastic regime, the scaled contact pressure p/G lies on one single curve, again underlying the 
importance of the shear modulus for the indentation behaviour. The only exception is given by the softest material. 

In this case the energy difference between the fee and the hep structures is rather small, and we observe instabilities 
at the relaxed free surface of the crystallite, which are reffected in the hardness curve. Fig. 5. 

C. Plasticity 

Figs, 6 and 7 give an atomistic presentation of the dislocations which developed in the substrate after the onset 
of plasticity. We analyzed the local atomic structures using a method based on the angular correlation of nearest- 
neighbour atoms. We found this method to be superior for our purposes than the more common analysis based on 
the centrosymmetry parameter. Only atoms deviating from the fee structure are visualized: the hep stacking faults 
(red) are surrounded by unidentifiable structures of low symmetry (grey) and a very small amount of particles in a 
bee coordination (green). These grey and green atoms occur at the boundaries of the dislocation loops due to the 
strong lattice deformations existing there and mark the ends of the loops. The free surface is also coloured grey. 

Fig. 6 displays the plasticity shortly after the nucleation of the first dislocations; note that the two materials shown 
here have a similar yield point (cf. Fig. 4), so that the amount of dislocated material can be compared. The damage 
is more concentrated around the indenting sphere for the stiffer material. Note also that the size of the dislocation 
loops has increased for the softer material. 

For fully developed plasticity (Fig. 7) we observe how prismatic dislocation loops have formed and were driven 
away from the indenter. In all cases which we studied, the (112) {111} glide system was activated first, and later, 
i.e. under larger stress gradients, also the (110) {111} glide system. This is a mere consequence of crystallography 
and is refiected by the smaller Schmid factor for the latter glide system; as expected, the changed elastic stiffness in 
our model crystals does not influence which glide systems are activated. However, it does affect the form and size of 
the plastic zones. Thus, for smaller elastic stiffness, (Fig. 7a), evidently more and smaller loops have been generated. 
This is in agreement with our finding that for smaller stiffness, also the unstable stacking fault energy - and hence the 
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barrier to dislocation formation and slip - is smaller: hence we have more loops. In contrast, for the stiffcr material 
(Fig. 7b), the nucleation of loops is retarded, and fewer but larger loops are formed. 



IV. CONCLUSIONS 



1. The Morse potential allows to systematically change the materials properties. We concentrate on changing the 
bulk modulus (and thus the elastic stiffness), while keeping the cohesive energy (bond strength) constant. In 
this potential, all elastic constants are changed in proportion to each other. 

2. Also the unstable stacking fault energy, which can be viewed as the resistance of a dislocation to gliding, changes 
in proportion to the elastic stiffness. The stable stacking fault energy, which is relevant for the width of partial 
dislocations, is quite unaffected from changes in the elastic stiffness. 

3. Not unexpectedly, the elastic part of the indentation curve scales with the clastic stiffness. 

4. The yield point, i.e., the indentation at which plasticity sets in, does not depend - apart from fluctuations - on 
the elastic stiffness. 

5. The load drop, which is experienced when plasticity sets in, increases in proportion to the elastic stiffness. 

6. The material hardness H is proportional to the theoretical strength, and thus on the elastic stiffness; it is 
ff= (0.20- 0.25) -G. 
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TABLE I: Fitted parameters of the potentials, D, a, ro, cf. Eq. (1). Materials properties determined from these potentials, B, 
G, 7„. 
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FIG. 1; Correlation between the shear modulus G and the bulk modulus B in the series of Morse potentials investigated. 
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FIG. 2: Generalized stacking fault energy as a function of displacement in [112] direction in the series of Morse potentials 
investigated. 



^ 700.0 

— 600.0 

500.0 

^ 400.0 

iS 300.0 

^ 200.0 

i£ 100.0 
^ 0.0 

^ -100.0 

I 0.0 100.0 200.0 300.0 400.0 500.0 600.0 

" Shear Modulus G (GPa) 

FIG. 3: Dependence of the stable and unstable stacking fault energies, 7„ and 7s on the shear modulus G in the series of Morse 
potentials investigated. 





FIG. 4: a) Dependence of the indentation force F on the indentation depth d in the series of Morse potentials investigated, b) 
Forces F normalized to the reduced elastic modulus Er- 




FIG. 5: Contact pressure p, normalized to the shear modulus G, as a function of indentation depth d in the series of Morse 
potentials investigated. After the onset of plasticity, the contact pressure defines the hardness of the material. 
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FIG. 6: Incipient plasticity determined at an indentation depth of ci = 9.7 A for an indentation target with bulk modulus of (a) 
B = 134.4 GPa h) B — 300.7 GPa. In each subfigure, the plots show a view from the bottom and three perspective side views, 
which can be identified by the axes drawn. The grey area shows the surface. The atoms colored in red display dislocation 
loops. 




(b) 



FIG. 7: Fully developed plasticity determined at an indentation depth of d = 11.9 A for the two cases shown in Fig. 6. 



